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Abstract The inelastic hard sphere model of granular 
material is simple, easily accessible to theory and simula- 
tion, and captures much of the physics of granular media. 
It has three drawbacks, all related to the approximation 
that collisions are instantaneous: 1) The number of colli- 
sions per unit time can diverge, i.e. the "inelastic collapse" 
can occur. 2) All interactions are binary, multiparticlc con- 
tacts cannot occur and 3) no static limit exists. We extend 
the inelastic hard sphere model by defining a duration of 
contact t c such that dissipation is allowed only if the time 
between contacts is larger than t c . We name this gener- 
alized model the TC model and discuss it using examples 
of dynamic and static systems. The contact duration used 
here does not change the instantaneous nature of the hard 
sphere contacts, but accounts for a reduced dissipation 
during "multiparticle contacts" . Kinetic and elastic ener- 
gies are defined as well as forces and stresses in the system. 
Finally, we present event-driven numerical simulations of 
situations far beyond the inelastic collapse, possible only 
with the TC model. 
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Introduction 

Granular media consist of discrete particles. Their interac- 
tion is governed by two major concepts: excluded volume 
and dissipation. Since the particles are solid, each particle 
occupies a certain amount of space, and no other parti- 
cle may enter this volume. If another particle approaches, 
the pair eventually collides. During collisions, energy is 
lost from those degrees of freedom (linear or rotational 
motion) which are important for the behavior of the ma- 
terial. Heat or sound are radiated and plastic deformation 
takes place so that energy is irreversibly lost Q. 

A model accounting for both excluded volume and 
dissipation is the inelastic hard-sphere (IHS) molecular 
dynamics with dissipative binary interactions. It is fre- 
quently used for the simulation of granular media, and 
is attractively simple. Between collisions, particles move 
freely through space. When two particles touch, their ve- 
locities are instantly replaced by new velocities calculated 
from a collision rule: 



w = c(u,n), 



(i) 



where U are the particles' velocities before the collision, 
and W are those after the collision. TZ denotes the parti- 
cles' positions at the time of collision. In theory, C could 
be anything, but in practice it is restricted by physical 
considerations, i.e. the particles may not interpenetrate, 
Galilean invariance, conservation of momentum, dissipa- 
tion of energy, etc.. U may also contain angular velocities, 
and C can be chosen to mimic real particles. Simulations 
using the hard sphere model can be very fast because to 
simulate a collision, the computer needs only to evaluate 
Eq. (|l|). On the other hand, if forces between particles 
are specified, the computer must integrate a differential 
equation over several time steps for each collision. Note 
that the assumption of an ideally hard potential is also 
used in kinetic theory and the Boltzmann or Enskog ap- 
proaches |2|-^) , which facilitates comparisons between sim- 
ulation and theory. All the collision rules we consider in 
this paper can be written 
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= Vi 2 ± 



[(v 2 - vi) -n] n, 



(2) 



where n is a unit vector joining the line of centers, and is 
the velocity of particle i. But the most important symbol 
appearing in Eq. (|^) is r, the restitution coefficient. The 
following equation for r can be derived from Eq. (0): 

-- ^- v 'i)' n (3 ) 



(v 2 - vi) • n 
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so r is the ratio of the component of the relative velocity 
along the line of centers after the collision to its value be- 
fore the collision. If r = 1, collisions conserve energy, and 
are said to be elastic. For < r < 1, energy is dissipated, 
and the collisions are inelastic. Usually r is considered to 
be a property of the material, and set to a constant which 
is the same for all collisions. 

The IHS model is best accessible to simulations and 
theory, but it has three general problems. First, there is 
the singularity of inelastic collapse: an infinite number of 
collisions can occur in finite time. Secondly, the collision 
rule treats only binary interactions, but in reality, many 
grains can interact, and these multiparticle interactions 
arc different from a sequence of binary collisions. Finally, 
there are no enduring contacts between particles, and no 
analog to various physical quantities, such as the energy 
stored in inter-particle contacts, exists. In this paper, we 
present the "TC model" , which is an extension of the IHS 
model that remedies these three problems. The collisions 
are still instantaneous, but we suppose that two parti- 
cles influence each other during a time t c after the col- 
lision. Specifically, if a particle experiences two collisions 
separated by a time less than t c , a multiparticle event is 
assumed to occur, and the second collision dissipates no 
energy. Except for this additional rule, the TC model is 
identical to the IHS model. 

The problems of the IHS model are discussed in more 
detail in Sec. II. In Sec. Ill, we review various attempts 
to solve these problems, including the TC model. Sec. IV 
describes the TC model in detail and applies it to the 
simple example of a particle lying on a flat surface. The 
elastic energy of a two-dimensional hard sphere gas is de- 
fined in section V and in section VI results on dissipative 
systems are presented that could be achieved by using the 
TC model whose consequences and future perspectives are 
discussed in section VII. 



Problems of the IHS model 

In this chapter, we discuss the problems of the IHS model 
that any improvement of it would have to correct. All 
three problems have essentially one origin: the potential 
used between the centers of mass of two colliding particles 
is unphysically stiff. The instantaneous collisions imply an 
interaction potential which is constant when there is no 
contact, and suddenly becomes infinite when the parti- 
cles touch. Therefore, momentum exchange takes place in 
zero time and thus the corresponding forces are infinite, 
however, acting for zero time only. In a real system the sit- 
uation is different: each contact takes a finite time during 
which large, but finite forces act. The infinitely stiff hard- 
sphere interaction is only an idealization or simplification 
of a smooth repulsive pair-potential. 

2.1 

Inelastic Collapse 

The most dramatic consequence of the infinitely stiff in- 
teraction potential used in the IHS model is inelastic col- 
lapse, which manifests itself as an infinite number of colli- 
sions in finite time. It was first discovered while studying 



the one-dimensional (ID) model system of a column of 
dissipative particles hitting a wall. The occurence of the 
inelastic collapse can be estimated using the product of 
the number of particles N and the dissipation per contact 
( 1 — r) . The effective dissipation £ = N(l — r) has a critical 
value of £ c w ir above which collapse occurs. The above 
value of £ c was calculated with the independent collision 
wave (ICW) model §. With slightly different arguments 
using the "cushion model" || , the value was evaluated as 
$ c wh[4/(l-r)]. The ICW model seems to work better in 
the inelastic limit, whereas the cushion model is superior 
in the elastic limit || . 

Inelastic collapse is also present in two dimensions (2D). 
In freely cooling systems a minimum of three particles is 
enough to lead to the collapse, if dissipation and density 
are large enough |],[l0|]. In larger assemblies, the inelas- 
tic collapse occurs, but it involves just a few particles 
arranged almost along a line. This leads to the conclu- 
sion that the inelastic collapse is mainly a ID effect [[ll] 
and that the one-dimensional predictions for the critical 
£ should work also in 2D. In fact, inelastic collapse in 2D 
unforced simulations can be predicted reasonably well by 
using the ID criteria with £ = (vl/d)(\ — r) (y is the 
fraction of the total area covered by the disks, I is the 
lenght of one side of the domain, and d is the disk diam- 
eter) |l^]. In a container in the presence of gravity, this 
expression for £ is equivalent to the number of layers of 
particles when the granular material is at rest. With these 
boundary conditions, the inelastic collapse likely occurs 
for small energy input and large £ [jl3l . Vibrated contain- 
ers with large filling heights cannot be simulated with the 
IHS model 

In two dimensions only a small fraction of the particles 
in the system is involved in the inelastic collapse. This im- 
plies that it is both physically insignificant for real particle 
assemblies and a major drawback for numerical simula- 
tions of dissipative systems. Therefore, any improvement 
of the IHS model must avoid the singularity of inelastic 
collapse. 



2.2 

Multiparticle interactions 

In the IHS model, true multiparticle interactions are im- 
possible because collisions are instantaneous. Multipar- 
ticle interaction is built up out of many two-particle in- 
teractions. In a real system the situation is different: Each 
contact takes a finite time so that multiparticle contacts 
are possible. The difference between two- and multiparti- 
cle contacts was examined for one- and two-dimensional 



model systems [15-171, and it was found that less energy 
is dissipated in multiparticle events than in an equivalent 
sequence of binary collisions. Any improvement of the IHS 
model should have the property that energy dissipation is 
reduced during multiparticle interactions. 



2.3 

No static limit 

Another problem with the IHS model is that the static 
limit does not exist, i.e. there is no way to represent endur- 
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ing contacts between particles. For example, in the frame- 
work of the IHS model, a particle cannot rest motionless 
on the ground. We disc uss this simple example in more 
detail in subsection 4.2, Another way to formulate this 



problem is by c ons idering the energies in the system, see 
also subsection 4.1, and especially the elastic contact en- 



ergy that is not defined in the IHS model. 

The translational kinetic energy E of the system is 
E = (1/2) Yli=i m i( v i + u l ) ; w ith the mass m.;, the fluc- 
tuation velocity Vi, and the flux velocity Ui of particle i. 
In the following we will mainly consider situations with 
Ui = 0. The potential energy E p is zero in absence of an 
external body force like e.g. gravity. In addition to E and 
E p , real materials have an elastic energy E e i at each con- 
tact, which is not defined in the context of the IHS model. 
In classical elastic systems, the total energy, i.e. the sum of 
kinetic, potential, and elastic energy, E + E p + E e i = E to t, 
is always conserved. In a dissipative system without a 
source of energy, the total energy tends towards a con- 
stant while the kinetic energy tends towards zero in the 
long time limit. This state is referred to as static, not to 
be confused with a "quasi-static" state, defi ned in the con- 



(h) 



text of the TC model below in subsection 4.1 



(hi) 



More specifically, consider a realistic granular material in- 
side a box, under the influence of gravity and in the ab- 
sence of energy sources. Due to dissipation, the material 
will loose energy. The potential and elastic energies will 
approach constants values while the kinetic energy tends 
towards zero. Eventually all particles are at rest, with 
E = 0, and the elastic energy is, as a rule, larger than 
zero, since a certain number of contacts is necessary to al- 
low for a stable static configuration. If the particles would 
be dissipative hard spheres, the inelastic collapse can oc- 
cur long before E vanishes so that the system will never 
reach a configuration with constant E to t- The inelastic 
collapse brings the system to an artificial halt. A static 
configuration with zero kinetic energy could be reached (iv). 
in a hard-sphere system by piling the spheres on top of 
each other (just in contact with no overlap). For this sit- 
uation E e i is not well defined and those touching hard 
spheres violate the rule of instantaneous contacts. Since 
in this artificial limit no elastic energy is defined, contact- 
forces and stresses are also not properly defined. This is 
an argument against the IHS model itself, which has by 
construction no static or "zero-temperature" limit. 
A method that loosens the restriction of instantaneous 
contacts and in return defines contact forces is the so- 
called contact dynamics, (see Ref. |l8) and references therein). 
The TC model, which we present in this paper, also pro- 
vides a way to define the elastic energy. 



Proposed modifications to the IHS model (v). 

In this section we review the various extensions and modi- 
fications of the IHS model which have been proposed. The 
main goal of these suggestions has been to remove inelas- 
tic collapse, since it is the most conspicuous problem of 
the IHS model. 

(i) . Particles with relative energy below a critical threshold 
can be merged into a "cluster" by setting their rela- 



tive velocity and separation to zero. If another parti- 
cle hits such a cluster, the momentum transport inside 
the cluster takes place instantaneously in the sequence 
of the largest relative velocity (LRV) fpTofj . With this 
method clusters can grow, and, given strong enough 
energy input, clusters may also be destroyed again. 
The deterministic LRV model has been successfully 
implemented in ID, but stresses and energies in the 
bulk are not defined. 

Several authors suggest a stochastic addition of trans- 
lational or rotational energy, as soon as the relative 
velocity after a collision drops below a critical value 
|L3],^(|. Also a small rotation of the relative velocity 
after contact of less than 5 degrees, seems to hinder 
the inelastic collapse [^(J, since correlations between 
successive collisions are diminished. Those stochastic 
models prevent the occurence of the inelastic collapse 
in the systems examined. However, neither is it clear 
whether the collapse can be circumvented under all 
conditions, nor has the physical relevance of the ran- 
dom energy input been discussed so far. 

Another method involves internal modes of every par- 
ticle. At each collision these modes may be agitated, 
and their energy is dissipated only on a time-scale 
longer than the duration of a contact. If a particle suf- 
fers an additional collision within this time, then en- 
ergy can be transferred from the internal modes back 
into translational motion. At least in a coolin g; sy stem 
in ID the inelastic collapse is prevented J2^,g2|, and 
possibly this model leads to an exp lana tion of the ran- 
dom energy input mentioned above | (ii)| . The drawback 
of this method is the effort necessary to model the in- 
ternal modes. 

A frequently used way to reduce dissipation in the low 
velocity regime (which has also been observed experi- 
mentally) is a velocity dependent restitution coefficient 
r(v ) based on the assumption of either viscoelastic 
jL5],^3| or plastic |^4],[2f| contacts. The velocity depen- 
dence seems to avoid the inelastic collapse in the ab- 
sence of walls and external forces [jl(|, however, its use 
for other boundary conditions resembling systems on 
earth, has not been examined up to now. Note that the 
dependence of r(v) on the collision velocity still con- 
cerns binary collisions with varying velocity, whereas 
the TC model discussed below concerns the transition 
from binary to multiparticle contacts - a completely 
different approach. Even when a velocity dependence 
of r can simply be added to the TC model, we avoid 
this in the following for the sake of simplicity. 

Instead of feeding energy into the system, another ap- 
proach just switches off dissipation following certain 
rules. The idea is to decide whether a particle still feels 
its previous collision partner when colliding with the 
next one. This is important since the presence of a 
third particle will affect the collision of a given pair. 
One can switch off dissipation if the next collision part- 
ner of a particle is detected within a critical distance 
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A c ]l2fl . In fact, it makes no sense to treat particles as 
separate objects if their surface-surface distance is of 
atomic size and, more technically, the numeric errors 
can become large when the distance between the par- 
ticles is many orders of magnitude smaller than the 
particle diameter. 

A similar argument for switching off dissipation is that 
two collisions cannot be treated as separate events if 
they take place within a too short time-interval t c that 
corresponds to the duration of a contact. It has been 
shown that collisions which overlap (in time) lead to 
weaker dissipation than a series of binary collisions 
(separated in time) |l5|,|l^,|6). Thus one sets the resti- 
tution coefficient to its elastic limit r = 1, if a par- 
ticle suffers more than one collision within time t c 
@|^,|§. With the time t$ , that passed by since 
the last collision n — 1 of particle i, the restitution 
coefficient for its n-th collision can be expressed as 



1 



for 
for 



> t c 
< t r 



(4) 



with < r < 1. Thus the type of a collision changes 
from inelastic to elastic when collisions occur too fre- 
quently. More specifically, a collision is elastic if at least 
one partner fulfills the above condition t n < t c . In 
general, r and t c can depend on the relative velocity 
and other parameters |l5|, p3[, so that the material's 
behavior can be adjusted app ropria tely using this de- 
pendence as discussed above | (iv) . Since Eq. (E) in- 
volves the duration of a contact t c we refer to this 
model as TC model in the following. 
Before turning to some more details of the TC model, 
we should mention that by using a velocity dependence 
for either the dissipation cut-off distance A c or the dis- 
sipation cut-off time t c , both can be connected: As- 
suming t c (v) oc v~ a , with the typical relative velocity 
v, leads to A c oc vt c (v) oc With a — 1 one has 

A c = const, and alternatively, with A c oc v one gets 
t c oc X c /v = const.. This unifies both the critical time 
and the critical length criteria into a single framework. 
In the following sections we will focus on the cut-off 
time (or contact duration) t c . 

To our knowledge, none of the models mentioned above 
has a sol id th eoretical background, except for the ap- 
proach in (iii) involving internal modes. The LRV method 
| (i)| has no reasonable static limit whe re stresses can be 
defined. The stochastic approaches (ii) require the choice 
of an a priori unknown sou rce o f fluctuation energy. The 
velocity dependence of r in (iv) was experimentally mea- 
sured for binary collisions, and is not necessarily impor- 
tant for multiparticle contacts. In the following we will 
discuss the physical relevance of the TC model, define its 
quasi-static limit, and finally apply it to selected exam- 
ples. 



The TC model in detail 

In Fig. [j], an interaction of two 'soft' particles (left) is 
compared to the interaction of two 'hard' particles (right). 



The words 'hard' and 'soft' correspond to non-smooth and 
smooth potentials between the centers of mass of two par- 
ticles respectively, and both are different approaches to 
model solid body interaction. Soft particles interact es- 
sentially during a time t c , indicated by the shaded region 
and the dashed vertical lines which mark the beginning 
and the ending of the contact. The contact duration t c 
of the soft particles is defined by these two instants of 
time. In the case of the hard particles, the interaction is 
instantaneous and the beginning and the ending coincide. 
However, in the TC model, the particles are considered to 
influence each other also during a time t c after each colli- 
sion (the shaded region). Note that the TC model in the 
limit t c = is identical to the IHS model and that both 
are identical to the elastic hard sphere model when r = 1. 



space 



space 




si 


















2' 


time 









Fig. 1. Schematic plot of the trajectories of two soft (left) and 
two hard (right) particles against time. The beginning and the 
ending of the interaction are marked by dashed and dotted ver- 
tical lines respectively and the time t c during which dissipation 
is affected is marked as shaded region. 



Assuming that t c is the physical duration of a contact, 
multiparticle contacts take place whenever a further colli- 
sion of the particles in Fig. l](b) occurs within the shaded 
area. This is the case when the typical time between in- 
stantaneous contacts t n gets smaller than the duration of 
a contact t c . The ratio r c = t n /t c is a measure for the 
existence of multiparticle contacts If r c ^> 1 pair in- 
teractions dominate, whereas for r c <C 1, one particle may 
be in contact with several others. Numerical simulations 
with various soft interaction potentials show that dissi- 
pation gets more and more ineffective with decreasing r c 



4.1 

The energies in the system 

In this subsection, we show how an elastic energy can be 
defined in the TC model. As was discussed above, there 
is no elastic energy in the IHS model. In a real system, 
only particles in contact with each other contribute to the 
elastic energy. As a consequence, we consider all particles 
which collided no longer than t c ago to be still in contact, 
and thus their energy contributes to a pseudo-elastic en- 
ergy E e . The translational kinetic energy E splits into a 
free kinetic energy Ek and the pseudo-elastic energy E e 
so that Ek = E — E e . Free means that the kinetic energy 
Ek can be dissipated whereas elastic energy E e cannot. In 
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both the IHS and TC models, E e i does not exist due to the 
non-smooth interaction potential, but in the TC model, it 
is replaced by E e . In table |l| we summarize the meaning 
of the symbols in the framework of the different models. 
From table |l| it is obvious that some type of elastic en- 
ergy is missing in the IHS model, while the TC model has 
three types of energy like a real system. However, a direct 
quantitative comparison of the energies in the framework 
of a realistic system to those in the TC model is far from 
the scope of this study |p9|] . 





Real 


IHS 


TC 


(free) kinetic energy 


E 


E 


Ek 


contact energy 


Eel 








elastic, kinetic energy 








E e 


potential energy 


E p 


E p 


Ep 



Table 1. Meaning of the symbols E, Ek, E e i, E e , and E v in 
the framework of a real system (or soft particle model), the 
classical hard sphere model (IHS) or the TC model. 



4.2 

The special case of one bouncing particle 

In order to discuss forces and stresses in the quasi-static 
limit, we examine the simplest possible case of a ball 
bouncing on a flat plane under the influence of the grav- 
itational acceleration g pointing in negative z direction. 
This system is essentially ID, and in the elastic limit, the 
particle will bounce forever. Introducing dissipation, with 
the restitution coefficient r, leads to a velocity just after 
the n-th collision v' n — —rv n as a function of the velocity 
v n just before. With the initial velocity v±, one has 

v n = r n -% . (5) 

The time between the collisions n and n + 1 is 

t n+ i = 2v n+ i/g , (6) 

for negative v„+i and negative g. In Fig. [2] the vertical 
position of a bouncing particle is plotted schematically. 
At each collision energy is lost and the particle stops at 
time t = t a . 



As already mentioned in subsection |2.3| the static limit 
of a real system is reached when E = 0, and E e i+E. p =const. 
The corresponding state of the TC model is, consequently, 
identified by Ek = 0, and E e + E p =const.. However, 
we denote it as the "quasi-static" limit, because E e > 
implies translational motion due to the kinetic nature of 
the elastic energy E e . A similar situation with E = 0, 
and E p =const. in the framework of the IHS model would 
mean that the system is artificially frozen. 

For example, consider a periodic system with realistic 
dissipative particles and without external forces (E p = 0) 
in its center of mass reference frame. Starting from an 
initial configuration with E > the system will evolve 
in time and the energy will decay. Note that this bound- 
ary conditon is t otal ly different from the situation dis- 
cussed in section |2~J| when gravity was active and walls 
were present. The only stable static situation (in the sense 
that E'tot =const.) would be one with all particles at rest 
(E = 0) and isolated, possibly just touching each other 
{E e i = 0). Such a situation one can denote as artificially 
frozen. Only in this case, kinetic and elastic energy vanish 
and the total energy can remain a constant. If two parti- 
cles would touch each other with E e i > 0, at least a part 
of their elastic energy will be transferred into relative ve- 
locity, eventually separating them, so that E > 0. This 
motion would eventually lead to more collisions reducing 
the total energy further. In other words: since no attrac- 
tive or external forces are involved, there is no reason for 
a stable overlap or deformation to exist in the long time 
limit. A similar argumentation for the TC model leads to 
the analogous conclusion Ek — E e = 0. In all other situ- 
ations the system evolves with time and the total energy 
is not constant. The IHS model, in contrast, comes to a 
halt when the inelastic collapse occurs and the long time 
limit cannot be reached. 



V 



n 



t 






— y — v — uk — > 



V 



n+l 



Fig. 2. Trajectory of a bouncing particle on a flat surface as 
a function of time. At time t s the particle is at rest. 



If the particle is infinitely rigid, as assumed by the 
IHS model, the particle will bounce an infinite number of 
times before t s . At times greater than t s , the IHS model 
is no longer defined. But this picture of an infinitely rigid 
bouncing ball makes no sense as soon as t n gets compara- 
ble to the duration of a contact t c . In that case the particle 
is in steady contact with the plate j30|. Therefore, in the 
TC model, r is set to 1 when t n < t c , and the particle 
bounces forever on the plate with a constant period which 
is less than t c . This is the TC model's representation of 
a particle lying on the plate. We now compare t s in the 
IHS and TC models. Starting with v x > gt c /(2r), the 
quasi-static limit is reached when t n+ \ < t c so that 



log r 



(7) 



In the hard-sphere limit n s — > 00 since t c — > 0. The time 
until the particle has lost all its kinetic energy is 



t (ms) 



^ tn+l 
n=l 



2vir 
9 



n=0 



1 



-ii 



(8) 



with ti from Eq. (j3|). The time until a particle with t c > 
reaches the quasi-static regime is 



t 



(TC) 



n=l 



t 



n+l 



r, n s -l 

~g~ ^ '' 

^ 71 = 



tf m \l-t c /h). (9) 



(i 



The difference At s between these two times is a measure 
for the difference between a soft particle (t c > 0) and a 
hard-particle (t c = 0) model: 



AU = t PHS) _ t (TC) 



1 



-tc 



(10) 



The last term in Eq. (|l(]) is obtained by assuming that the 
time between elastic collisions n > n s is approximately t c 
(what is almost true for r ~ 1). Note that At s will be 
small, because the contact time is usually small. 



4.2.1 

Performing averages 

In the framework of the TC model, an observable A has 
to be defined in average over a time-interval At: 



A{t) = (A) = (A) At = 



1 

At 



A{t') At'. 



(11) 



t-At 



Alternatively, ensemble averages can be performed, how- 
ever, this option will not be discussed here. The average 
makes sense only if it averages at least over the duration 
of a contact t c , since the TC model simplifies the real- 
ity during times smaller than t c . Therefore, averages over 
longer intervals should be taken to level out the details of 
the basic assumptions introduced in e.g. Eq. (E). 

We now discuss this average, using as an example one 
particle resting on a flat surface. A real, soft particle rest- 
ing on the bottom is represented in the TC model as an 
elastic, hard particle bouncing on the surface with a pe- 
riod t n < t c . Since it performs a periodic orbit with du- 
ration t n , one can set At — t n . Thus, integration is per- 
formed over one parabola of free flight of the particle. The 
mean velocity of the bouncing particle is u = (v) = 0, 
as expected for the quasi-static limit, the mean squared 
fluctuation velocity is {(v - it) 2 ) = (v 2 ) — {\/\2){gt n ) 2 , 
and the mean separation of the particle from the bot- 
tom is (z) — {\/\2)gt 2 l . With these quantities we may 
identify the energies defined above. The potential energy 
is connected to the separation from the bottom E p = 
nig (z) = (1 / '12) m(gt n ) 2 , disregarding an additive con- 
stant. The total translational energy is E = Ek + E e , with 
E = (m/2)(v 2 ) = (l/24)m{gt n ) 2 . Now, we identify the 
elastic energy with the kinetic energy of the particle(s) 
which suffered a collision no longer than t c ago. Since the 
particle collides with a rate i" 1 > t~ l , all its kinetic en- 
ergy contributes to E e so that Ek = 0. Note that the 
values Ek — 0, E p > 0, and E e > correspond to the 
quasi-static limit discussed above. 

In addition, one can calculate the force which the par- 
ticle exerts onto the bottom as the momentum exchange 
per unit time / = (Ap) — mg, as to be expected for a 
particle with mass m in the gravitational field p7[] . 



4.2.2 

The link to a linear elastic particle 

A soft, elastic particle in contact with the bottom has - in 
the framework of the simplest linear model |3l[] - the elas- 
tic energy V(S) = (l/2)fc<5 2 , with stiffness fc, and overlap 
or deformation S. At rest, it exerts the force fk = kSo = 



mg onto the bottom, so that the overlap or deformation 
is do = mg/k. When bouncing, its contact duration is 

tf = Tl/U0 = TTy/m/k , (12) 

what leads to the identity 5q = g{tf / tt) 2 . The elastic en- 
ergy of the particle at rest is thus V(So) — m(gtf) 2 /{2'K 2 ). 

Comparing the soft particle with the TC model from 
the previous subsection by using either of the relations 
(z) = So or E e = V(So), leads to 

t n = VUtf/n . (13) 

From the beginning of this section we remember that the 
particle reaches its quasi-static limit when t n < t c . Thus 
the contact duration t c w t n can be identified with the 
contact duration of the linear soft-sphere model tf , when 



disregarding the constant factor \JY2/t{ m 1. 
4.2.3 

A Gedanken-Experiment 

In order to clarify the meaning of the different energies 
calculated above, we assume that we are able to switch off 
gravity at any time during the period t n . A real particle 
lying on a table with zero kinetic energy will then begin 
to rise due to the elastic energy stored in the contact. 
Ideally, for r = 1, all elastic energy will be transferred 
to translational motion. Within the framework of the TC 
model, as discussed here, one can calculate the velocity a 
particle will eventually reach, after g is set to zero. Since 
the particle- velocity is phase-dependent, one has to per- 
form the average over all possible phases at which gravity 
might be switched off. Thus the integration over the pe- 
riod t n has to be split: During the first half-period for 
t < t n /2, the particle will keep its upwards velocity and 
move away from the bottom. In the second half-period 
the particle will suffer one collision with the bottom be- 
fore it moves upwards. Performing the integration one gets 
( v g~>o) = (l/ 24 )(s*n) 2 (l + r )- Note that the velocity af- 
ter switching off gravity, u g ^o , depends on the time when 
gravity is switched off. This reflects the fact that a too fine 
resolution in time uncovers some details of the simplifica- 
tions in the TC model. However, the elastic case r -> 1 
corresponds to the case when all elastic energy E e from 
the quasi-static regime is transferred into kinetic energy 
Ek- 



The TC model in elastic systems 

In the following we will mainly focus on the elastic limit 
r — > 1, and define the properties of interest like the stress 
tensor, the equation of state, the collision rate, and the 
elastic energy. For the simulations in this section we use 
an event driven (ED) simulation method as introduced by 
Lubachevsky |32|. The system has periodic boundaries, 



and neither walls nor gravity are present. The following 
discussion concerns only systems in equilibrium and in 
their center of mass reference frame. Strictly speaking, all 
collisions in the simulations discussed in this section are 
"elastic", because r = 1. However, the TC model distin- 
guishes between those collisions which indicate multipar- 
ticle events (where one of the partners has had a collision 
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no longer than t c ago), and the rest of the collisions which 
are truly binary. In this section, only the former type of 
collisions are called "elastic" , for they would also be elas- 
tic (in the TC model) when r ^ 1. We choose r = 1 here 
because this gives systems in equilibrium and allows to 
obtain good statistics by taking long time averages. 



5.1 

The stress tensor 

The stress tensor, defined for a test-volume V c , can be 
written component-wise as 



&af3 



1 

Vr. 



(14) 



The first sum runs over all contact points j, and the sec- 
ond sum runs over all particles i, both within V c [|33|,|34|]. 
The indices a and j3 denote the Cartesian coordinates, l a 
are the components of the vector from the center of mass 
of a particle to its contact point j, where a force with com- 
ponents acts. A particle i has a mass m, and a velocity 
with the components v a . 

In the static limit, the second term drops out, since 
all velocities vanish. In a dilute system without perma- 
nent contacts, the first term would be negligible. On the 
other hand, for a hard-sphere gas, the first term has to 
be treated differently, since no forces are defined. The dy- 
namic equivalent to fn is the change of momentum per 
unit time App/At Q. For a hard-sphere gas the stress 
due to collisions may be evaluated as an average over all 
collisions in the time interval At 



v a p{t) = {o a p) At = 

—j2 e MA Pf3 (t n )-j2 



(15) 



Here, the first sum runs over all collisions n occuring in 
the time between t — At and t. In general, the volume V c 
and the time-interval At have to be chosen large enough 
to allow averages over enough particles and enough col- 
lisions, but also small enough to resolve inhomogeneities 
in space and variations in time. Note that the result may 
depend on the choice of At and V c ||, so that this aver- 
aging procedure is not necessarily the best choice under 
all circumstances. 



5.2 

The equation of state 

In order to test the averaged stresses defined above, we 
compute the mean pressure P = (<7i + 02)12 from the 
eigenvalues o\ and oi of the stress tensor. The results of 
simulations with elastic particles, r = 1, and different vol- 
ume fractions v are compared to the pressure obtained by 
the two-dimensional equivalent of the Carnahan-Starling 
formula ||,^!|. In Fig. ||, we plot the reduced, dimension- 
less pressure 



PV 

P - 1 = — - 1 = 2vg{2a) 



against the volume fraction v — Nna 2 /V pi| ]. We use the 
pressure P, the volume of the system V, the total trans- 
lational energy E = (l/2)U^ 1 m,ivf, and the particle- 
particle correlation function evaluated at contact 



g(2a) 



(17) 



taken from Eqs. (28) and (110) in Ref. g. The agreement 
between theory and simulations is perfect since r = 1, 
already for a tiny system with N = 42. The simulations 
deviate from theory only for volume fractions larger than 
about 0.7 because then a solid state exists, i.e. the motion 
of the particles is hindered by their neighbors j3fj|j . The 
smallest system deviates slightly from the larger ones, in- 
dicating finite size effects in the crystalline, high density 
regime. However, we will not discuss the transition from 
a fluid, disordered to a 'solid', ordered system here. 
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Fig. 3. Reduced dimensionless pressure Po — 1 plotted against 
volume fraction v. Theory (solid line) is compared to simula- 
tions (circles) from a system with N = 1435 particles. In the 
inset simulations of a smaller and a larger system, with N = 42 
(small squares) and iV = 16680 (small triangles), respectively, 
are presented in addition. At least three simulations were per- 
formed for each v value and plotted over each other so that 
smeared out symbols indicate comparatively large fluctuations 
at v k, 0.7. 



5.3 

The collision rate 

The next quantity of interest is the collision rate, i.e. the 
number of collisions per particle per unit time. Thus, we 
define the collision rate in the simulations as the inverse 
of the typical time between contacts 



(16) c r = t; n 



2C 
N At 



(18) 
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with C = (C), the number of collisions per averaging time 
At. Note that we use C as the number of collisions within 
a time At, whereas the total number of collisions since 
the beginning of the simulation is denoted as Ct later on. 
The prefactor '2' stems from the fact that each collision 
involves two particles. In Fig. ^ we compare the simula- 
tion results with the Enskog collision frequency in 2D, as 
expressed in Ref. M, 



= 4a- 



N 



E , n s 
7T— g{2a) 



v\ m - -o(Po-i), (19) 

with the particle radius a, the total fluctuation kinetic en- 
ergy E, the total mass M, and g(2a) as defined in the 
previous subsection. From Fig. ^ we observe again a per- 
fect agreement between theory and simulation (circles) for 
v < 0.7. The difference between Eqs. (pGj) and ( JlS| ) is the 
constant prefactor oj = \/^7^( v t / a ) , with the thermal 
fluctuation velocity vt = y/2E/M. Note that both Ej- 
and E e contribute to E and thus to the collision rate. 
Therefore, C r can be split into a dynamic and a quasi- 
static contribution. The latter, the collision rate of elastic 
collisions is displayed in Fig. [| in addition to the over- 
all collision rate. In an elastic collision, at least one col- 
lisionpartner had a collision no longer than t c ago, see 
Eq. (||). Therefore, smaller t c values lead to lower collision 
rates (the decadic logarithm of t c = 10~ 3 s, 10 ~ 4 s, 10 ~ 5 s, 
10 _6 s, and 10 _7 s is given in the inset). 



10000 



I 1000 

c 

o 



o 
o 




Fig. 4. Collision rates plotted against volume fraction v from 
theory (solid line) and the simulations (circles) from Fig. |^. 
The dashed lines are inserted between the measured elastic 
collision rates (small symbols) to guide the eye. The system 
dependent prefactor of these simulations is u>o = 656.03 s _1 . 
The small symbols indicate the collision rates of elastic col- 
lisions according to the TC model with log 10 t c given in the 
inset. 



5.4 

Elastic particles, collisions, and energies 

In this subsection, we are interested in n e = N e /N, the 
fraction of particles that are elastic, i.e. which had a colli- 
sion no longer than t c ago. Furthermore, we want to esti- 
mate c e = C e /C, the fraction of collisions in which elastic 
particles participate, related to the elastic collision rate 
in the pre viou s subsection. Finally, following the ideas in 
subsection hi , we will split the total translational energy 
into a kinetic and an elastic part, i.e. E = Ek + E e and 
determine e e = E e /E, the fraction of elastic energy in the 
system. 

First we estimate the probability p(t c ) = 1 — N e /N 
that a particle had no collision since a time t c |37]]. We 
know that p(0) — 1, since no particle can suffer a collision 
in zero time, and p(oo) — 0, since any particle will even- 
tually collide (given that a > and that the system is not 
artificially frozen). Furthermore, we also know the prob- 
ability dt that a particle will collide within time dt. 
Thus we have the probability 1 — t^dt that the particle 
will not collide. Multiplying the probability that it did not 
collide until t c with 1 — t^dt finally gives the probability 
that it does not collide until t c + dt. First order Taylor 
expansion of p(t c + dt) — p(t c )(l — t~^}dt) around t c leads 
to 4zp(t) = pit)^ 1 . Integration for constant collision fre- 
quency t~j^ gives the fraction of 'elastic' particles 

n e = l-p(t c ) = l-cxp{-t c t E 1 ) . (20) 

In Fig. H we present the quality factor q n = (N e ) / (Nn e ) , 
the ratio of the measured fraction of elastic particles and 
of the analytical expression from Eq. (^). Each point cor- 
responds to an average over 30 snapshots from simulations 
with N = 1435 particles. 




Fig. 5. The quality factor q n , i.e. the ratio of (N e ) /N and 
n e from Eq. (pOK plotted against v. The simulations are the 
same as in Fig. S. Values close to unity indicate good agree- 
ment between theory and simulations for the t c values used. 
To distinguish between different t c values, log 10 t c is given in 
inset. 
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Data for different t c values can be obtained from the 
same simulation since each particle 'remembers' when it 
had its last collision at time tn and because no dissipation 
is active. Dissipation would make the sequence of collisions 
dependent on t c . The agreement between simulations and 
the theoretical prediction, Eq. (pfj|), is reasonable (q„ w 1). 
Deviations occur for small volume fractions v due to bad 
statistics and for large volume fractions v > 0.7 due to 
the transition to the ordered system. 

However, Eq. (^) is not the fraction of elastic colli- 
sions c e which instead has to be computed as the sum of 
probabilities that either both or only one of the collision 
partners belongs to the elastic particles. Thus one gets 



= n\ + 2n e (l - 7i e ) = 1 - exp(-2i c t/) , 



(21) 



almost in agreement with the simulations as displayed 
in Fig. [|. The quality factor q c = (C e ) / ((C) c e ) plotted 
against v indicates a difference between the measured and 
the calculated fraction of elastic collisions. Interestingly, 
the simulation data show that there occur about five per- 
cent less elastic collisions than expected from Eq. (gjj). 
The data presented in Fig. || were obtained by averaging 
(C e ) / (C) over a time interval At large enough to allow 
for more than about 10 4 collisions per particle. A more 
rigorus theoretical treatment that would lead to the exact 
value of c e is beyond the scope of this study and will be 
discussed elsewhere [B8|. 



the numerical simulations can, however, be understood 
via qualitative arguments. Particles with greater veloci- 
ties have a higher probability to collide and, therefore, 
have a greater collision rate than slower particles. Due to 
the greater collision rate, fast particles are more likely to 
contribute to E e and E e is increased since faster particles 
contribute with a greater energy Q . 
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Fig. 6. The ratio of simulational and theoretical results on 
the probability for elastic collisions plotted against v. Exact 
agreement would correspond to q c = 1. 



Finally, we estimate the elastic energy contained in 
the system in a simple minded, mean-field way. The frac- 
tion of elastic energy e e should be the product of n e and 
E when all particles would have the same (mean) energy 
and would contribute to E e with the same probability. Un- 
fortunately, the simulation results in Fig. ^ indicate that 
the mean-field approach is not valid. The discrepancy be- 
tween the mean-field estimate presented above for n e and 



Fig. 7. The ratio of simulational and theoretical results on 
the fraction of elastic energy in the system plotted against v. 
A quality factor q e — 1 would correspond to exact agreement. 
The value 5/4 was taken from a more elaborate kinetic theory 
calculation IBal. 



The TC model in dissipative systems 

One of the phenomena which received a great deal of at- 
tention in the last years is the clustering instability in 
rather dilute systems of dissipative particles p2],|39|,^0[. 
Initially given a homogeneous density and a Maxwellian 
velocity-distribution, the system cools due to dissipation 
[|l2],^(|. This cooling regime is unstable to perturbations 
with large enough wavelength, i.e. small enough wavenum- 
ber The homogeneous state can be well described 

by hydrodynamic theory, but the standard description 
breaks down as soon as the perturbations grow - assump- 
tions like homogeneity or "molecular chaos" are not longer 
true fio). Furthermore, it has been recognized that a 
hydrodynamic description breaks down due to the diver- 
gence of the collision rate and the connected dramatic de- 
crease of free volume during the inelastic collapse (i^, [l3| . 
The instability may be understood in a qualitative man- 
ner: The homogeneous state contains thermal velocity fluc- 
tuations. A convergent velocity fluctuation leads to in- 
creasing densities in certain regions. As the collision rate 
increases in these regions, so does the energy dissipation 
rate. If the energy dissipation rate is great enough, the 
pressure cannot reverse the convergent flow. If this process 
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is not terminated by sufficiently strong perturbations, it 
is self-stabilizing, causes clusters, and may eventually lead 
to the inelastic collapse. 

The clustering instability and the inelastic collapse were 
carefully examined in ID [@,|,|J,[0|-f47| and in 2D §, 
|l|,|^[^,|4§-[5l). Cluster growth could be described the- 
oretically in the case of irreversible aggregation (r = 0) 
Jf9|,|5^] and, more recently also a theory for the growth 
of density fluctuations was proposed for r ^= |53j. De- 
tailed examination of the inelastic collapse by McNamara 
and Young [|12j led to the picture of different 'phases'. In 
a periodic system without external forcing exists a criti- 
cal dissipation - connected to system size, volume fraction 
and restitution coefficient - above which clustering occurs 
and below which the system stays in molecular chaos. In 
the transient regime shearing modes or large scale eddies 
are frequently observed. The case of homogeneous cooling 
is rather well understood 28 1 so that we mainly focus on 
situations with rather strong dissipation when the system 
is no longer homogeneous. 

Using event-driven simulations, periodic 2D systems 
with side-length L = I /(2a) in 2D are examined in the 
following. A system contains N particles of radius a and 
volume fraction v = NTr(a/l) 2 = (tt / '4)iV '/ L 2 . In this sec- 
tion the particles are dissipative with a restitution coeffi- 
cient r. We apply the TC model as described above, i.e. we 
use Eq. (^) with the parameter t c to be specified. Initially 
we arrange the particles on an ordered lattice and give 
each of them a random velocity; then the system is equili- 
brated with r = 1 until a Maxwellian velocity distribution 
is obtained. Finally, dissipation is switched on and the 
simulation starts at t = s. 



6.1 

Inhomogeneous cooling - Finite size effects 

In this subsection we discuss typical simulations with vol- 
ume fraction v = 0.227, restitution coefficient r = 0.40, 
and contact duration t c = 10~ 5 s. The system size is var- 
ied so that N = 97776, 22960, 5740, 1435, 378, and 42 
particles fit into the system. In Fig. H, the dimension- 
less kinetic energy T = E(t)/E(0) is plotted against the 
rescaled time r = (0)t, with the initial collision rate 
t~ 1 (0) w 444 s" 1 for all simulations. 

Following the estimate of a critical restitution coeffi- 
cient, see Eq. (9) in Ref. [Q, below which the inelastic 
collapse can be expected, we compute for our simulations 



r c (N) « tan 2 



1 - 



1 



A 



pt 



(22) 



with the non-dimensional optical depth A op t = V nNv/2. 
Comparing the critical restitution coefficients r c (97776) = 
0.976, r c (22960) = 0.952, r c (5740) = 0.906, r c (1435) = 
0.821, r c (378) = 0.680 and r c (42) = 0.296 to r = 0.40 used 
for the simulations, we note that all are larger, r c > r, ex- 
cept for the tiny system with N = 42. However, due to the 
TC model used in our simulation no collapse is observed. 
The dashed line in Fig. || gives the analytical solution for 
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Fig. 8. T as function of r for simulations with N = 97776, 
22960, 5740, 1435, 378, and 42, v = 0.227, r = 0.4, and t c = 
10 s. The dashed line gives the solution of the homogeneous 
cooling state Eq. (E3). 
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Fig. 9. Collision rate C r , see Eq. (|l8|), as function of r from 
the same simulations as in Fig. H. The dashed line gives the 
solution of the homogeneous cooling state. 



the homogeneous cooling state (HCS) of smooth particles 
(see H and references therein): 



T= 1 



1 - r A 



(23) 



Only the smallest system is close to the theoretical pre- 
diction, all others deviate stronger with increasing system 
size. 

The collision rate in the homogeneous cooling state is 
directly linked to the fluctuation velocity vt, see Eq. [l9| 
Because vt is the only quantity that changes during the 
simulations, we can also express the collision rate in terms 
of \/T = vt/vq so that due to normalization Tq = 1, one 
has C r = C r (0)VT. In order to test this prediction, we plot 
C r against r in Fig. M from the same simulations as in Fig. 
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Fig. 11. Snapshots at time t — 5 s, i.e. r = 2222, from ED simulations started from identical initial conditions, with N = 5740, 
L — 140, v — 0.227, r — 0.4, and different contact duration t c . The mean fluctuation velocity (relative to the center of mass) 
is color-coded red (vt/vo > 0.045), green (vt/vq ~ 0.02), and blue [vt/vq < 0.01), with initial mean velocity vo = vt(0) = 
0.3615 ms -1 . 



As can be seen in both figures land§ for small r < 5 
the simulations agree with the theory. For larger r de- 
viations from the homogeneous cooling regime occur and 
the energy decay slows down due to the density instability 
and the build-up of clusters. The behavior of the energy 
as a function of time is independent of the system size 
up to t w 200 when the small system TV = 378 starts to 
deviate from the larger systems. The deviation from the 



common behavior occurs later with increasing system size, 
indicating finite-size effects. A closer examination of snap- 
shots from the simulations leads to the conclusion that the 
deviation from the slow cooling regime occurs when the 
clusters have reached the system size. The tiny system 
N — 42 is anyway too small to allow large scale structures 
and therefore follows the HCS prediction more closely. 
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In Fig. JlG snapshots of the simulation with N = 22960 
from Fig. Hare displayed. With increasing time t struc- 
tures build up in the system and grow in size. In the red 
and green regions in the centers of the clusters the col- 
lision rate is largest. The size of the clusters grows with 
time and the qualitative cooling behavior of the system 
changes as soon as the clusters get to be as large as the 
system. 



6.2 

Variation of the contact duration 

In the following set of simulations, t c is varied while the 
other parameters are fixed to N = 5740, v w 0.227, and 
r = 0.4. In Fig. [ll] snapshots from the simulations with 
different contact duration are displayed. The simulations 
were set up with identical initial conditions and the snap- 
shots were taken at the same time. For t c = 10 _2 s cool- 
ing is almost hindered, since the time between collisions 
is much greater than the dissipation cut-off t c . The simu- 
lations with t c < 10 -3 s all show density instabilities and 
clusters, however, the pictures differ in details. The change 
in behavior between t c = 10~ 2 s and 10 -3 s will be di- 
cussed in more detail elsewhere [fj8[ . Only for very small 
t c < 10~ 8 s the simulations lead to almost identical re- 
sults. The color of the snapshots represents the particle 
velocities in the center of mass reference frame. A large 
cluster with different colors is thus not a block without 
internal motion but rather a liquid-like arrangement with 
strong internal shearing. 

In Fig. [l2|, the dimensionless kinetic energy T is plot- 
ted against the rescaled time r = (0) t, with the initial 
collision rate i^ 1 (0) ~ 444 s _1 . As contact duration, var- 
ious values between t c = 10~ 2 s, and 10 _12 s were used 
(log 10 t c is given in the inset for identification). The sim- 
ulations with t c < 5 x 10~ 4 s lead to the same functional 
behavior of T as function of r, except for small deviations 
for large r, as can be obtained from the inset. However, 
the data with t c < 10 -6 s cannot be distinguished, i.e. for 
small enough t c the energy of the system is not influenced 
by the contact duration. For very small r < 5, the simu- 
lations agree with the theoretical result for the homoge- 
neous cooling state from Eq. (p3|) . For larger r deviations 
from the homogeneous cooling regime occur and the en- 
ergy decay slows down due to the density instability and 
the build-up of clusters. 

The simulations with the largest t c cool slowest since 
the number of elastic collisions increases with tc, see Eq. 
(H|). Inserting t c = 10~ 3 s and t^}(0) into Eq. @, which 
was derived for elastic system in equilibrium however, 
leads to c e « 0.59, whereas the contact duration of t c = 
10~ 4 s leads to a much smaller fraction of elastic con- 
tatcs c e w 0.085. Evidently, dissipation is almost inactive 
when t c is as large as t c — 10~ 2 s so that c e w 0.9999. 
On the other hand, the effect of t c will be negligible for 
t c < 10~ 5 s when the fraction of elastic collisions vanishes 
c e < 9 x 10~ 4 . In summary, the contact duration t c slows 
down dissipation in the system - the larger t c the stronger 
the effect. For small enough t c the system is not affected 
and the variation of the energy T with the contact du- 
ration t c is rather weak. Note that other quantities, as 
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Fig. 12. Dimensionless energy T plotted against r from simu- 
lations with r — 0.4 and log 10 t c as given in the plot. The inset 
is a zoom into the plot at large r and small T. The dashed line 
corresponds to Eq. (123) with r = 0.4. 
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Fig. 13. T plotted against Ct/N for the same simulations as 
in Fig. |l2[ 



e.g. the total number of collisions per particle Ct/N, can 
vary strongly with t c , as evidenced in Fig. [l3| The smaller 
t c the more collisions occur in the system during a fixed 
time interval. Sometimes, for t c < 10~ 10 s, jumps in Ct/N 
are observed. As a consequence, the total number of col- 
lisions looses its meaning as a system-inherent time scale 
as soon as the system becomes inhomogeneous. This can 
also be understood when recalling that an inhomogeneous 
system is a system that has non-constant density, fluctua- 
tion velocity, pressure, and collision rate. Different parts of 
the system, evolving with different collision rates, cannot 
be assumed to follow a common time-scale. On the other 
hand, the variation of Ct/N with t c allows to see the TC 
model also as a way to reduce the computational effort 
by decreasing the global number of events that have to be 
handled, however, without affecting physical observables 
like T. 




Fig. 14. (a) The fraction of elastic collisions C e /C as function of r with r = 0.4 and log 10 t c as given in the inset. The simulations 
are the same as in Fig. |l^, however, each point represents an average over all collisions in the time interval At = 0.1 s. (b) The 
fraction of elastic particles N e /N plotted against r for the same simulations as in (a). Each data point represents an average 
over 20 snapshots in the time interval At = 0.1s. Note the different vertical axis scaling in (a) and (b). (c) The fraction of 
elastic energy E e /E plotted against r for the same simulations and the same averaging procedure as in (b). 



Another way to report the effect of the TC model on the 
simulation results is to take a closer look at the fractions 
of elastic collisions, particles, and energy in the system. 
In Fig. |lj(a), (b), and (c) these quantities are displayed 
respectively. The fraction of elastic collisions c e = C e /C 
is systematically larger than both the fraction of elastic 
particles n e = N e /N and the fraction of elastic energy in 
the system e e = E e /E. The latter two quantities are al- 
ways comparable, whereas c e depends on t c in a different 
way. The fraction of elastic collisions is correlated to t c , 
at the beginning of the simulation, as long as the system 
is homogeneous. For larger r the fraction of elastic colli- 
sions fluctuates around at a finite value, independent how 
small t c is. Even the simulation with t c = 10 _12 s has a 
rather large fraction of elastic collisions. However, setting 
a large C e in relation to the large total number of collisions 
Ct/N, or cquivalcntly C r , leads to the conclusion that the 
TC model adjusts c e to an finite, almost constant value. 
In contrast, the elastic energy and the fraction of elastic 
particles decrease systematically with t c . These two quan- 
tities are directly affected by t c . Concerning averages we 
must remark that the procedure to obtain C is different 
from the one used to obtain average particle numbers or 
energies. While all collisions in the averaging time interval 
are summed up to C, only a certain number of snapshots 
is evaluated to compute (N e ) and (E e ). 

6.3 

Variation of the restitution coefficient 

For the parameter study in this subsection, the system 
with N = 5740, v = 0.227, and a fixed contact duration 
t c = 10^ 5 s is used. This choice of t c is arbitrary, however, 
as we have shown in the previous subsection, it is a rea- 
sonable choice in the framework of the TC model: t c has 
a negligible effect on T. Furthermore, such a t c value is of 
the same order of magnitude as the contact duration of 
e.g. two steel spheres with radius a = 1.5mm fjlj. 



LOO 




10 100 1000 

X 



Fig. 15. Dimensionless energy T plotted against r from sim- 
ulations with t c — 10~ 5 s and r as given in the plot. The inset 
is a zoom into the plot at large r and small T. The points are 
simulation data, the lines give Eq. (ji^) for the corresponding 
r values. 



In Fig. |15|, the dimensionless kinetic energy T is again 
plotted against the rescaled time r, with the initial col- 
lision rate t^}(0) s=s 444 s -1 . Simulations are performed 
for different restitution coefficients as given in the plot, 
and are compared to the analytical solution of the homo- 
geneous cooling state, see Eq. (p3|). Only for r = 0.99, 
we evidence reasonable agreement with the theory. With 
decreasing r, deviations from the theoretical curve occur 
already for smaller and smaller r. However, for r < 0.4, 
there are only small differences between simulations with 
different r, as can be seen more clearly in the inset. 
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Fig. 16. (a) The fraction of elastic collisions C e /C as function of r with t c — 10 -5 s and r as given in the inset. The simulations 
are the same as in Fig. |l5|, however, each point represents an average over all collisions in the time interval At — 0.2 s. (b) The 
fraction of elastic particles N e /N plotted against r for the same simulations as in (a). Each data point represents an average 
over 20 snapshots in the time interval At = 0.2 s. Note the different vertical axis scaling in (a) and (b). (c) The fraction of 
elastic energy E e /E plotted against r for the same simulations and the same averaging procedure as in (b). 



We do not present a plot of T against C t /N here, 
but make some qualitative remarks on the total num- 
ber of collisions. We evidence that the total number of 
collisions per particle varies with r as it varies with t c . 
For large r and small Ct/N the energy behaves as T = 
exp[— 5(1 — r 2 )Ct/N], as can be simply derived using Eq. 
( f23| ) and integrating the collision rate C r over r. This be- 
havior is obtained only for r « 1, already for r < 0.95 the 
simulations deviate from the theoretical prediction. 

For strong dissipation, only a few particles perform 
many collisions. One can expect that the time between 
the collisions of such particles drops below the threshold 
t c so that r is set to unity according to Eq. (Q). The TC 
model is active, hinders dissipation, and thus evades the 
inelastic collapse. The number of particles affected by the 
TC model will be discussed in the following. In Fig. [l6] 
we present data on c e , n e , and e e in a way similar to Fig. 
|l4| . Again, we obtain that c e behaves differently from n e 
and e e . All simulations with r > 0.80 have a negligible 
fraction of elastic collisions (c e < 2 x 10~ 3 for long times). 
The simulations with stronger dissipation (r < 0.80) have 
an increasing number of elastic collisions with increasing 
dissipation, i.e. decreasing r. 

From the data on N e , we can estimate the number of 
particles with the largest collision rate, which typically 
are involved into elastic collisons. A fraction of n e = 0.002 
corresponds in the case of N — 5740 to N e ss 10 particles. 
Even in the case of very strong dissipation only about 10 
particles are affected by the TC model. 

7 

Conclusion and Outlook 

In this study we discussed the TC model, an extension of 
the frequently used inelastic hard-sphere model. Introduc- 
ing the contact duration t c as a material parameter, multi- 
particle interactions are defined in some sense. They con- 
cern particles with large collision rates which are assumed 



to contribute to the elastic energy in the system (which 
cannot be dissipated). The TC model can reach a quasi- 
static situation when only elastic and potential energy 
remain. Dissipation is locally inactive for large collision 
rates, i.e. the elastic limit where multiparticle contacts 
occur, and active for rare events, i.e. the dissipative limit 
where contacts are binary almost always. The TC model 
allows simulations in ranges of parameter space, where the 
classical inelastic hard-sphere model breaks down due to 
the inelastic collapse. Potential, kinetic, and elastic ener- 
gies as well as stresses and forces are defined as averages 
over time-intervals comparable to t c . The material param- 
eter t c can be identified with the contact duration tf of a 
simple linear particle model, involving e.g. particle mass 
and stiffness. 

Furthermore, mean-field estimates for the fraction of 
elastic particles n e , the fraction of elastic collisions c e , 
and the fraction of elastic energy e e in the system are pre- 
sented. The simulation results indicate that a more elab- 
orate theory is required to explain the obtained discrep- 
ancies. However, the definition of n e , c e , and e e is also 
valid in non-equilibrium, dissipative systems. Detailed ex- 
aminations of the inhomogeneously cooling situation leads 
to the conclusion that the TC model affects a small frac- 
tion of the particles only - just as the inelastic collapse. 
Therefore, we beleive that the TC model removes inelastic 
collapse in a physically reasonable way, without disturb- 
ing the global behavior of the system, i.e. the clustering. 
This is strictly true for realistically small t c values, be- 
cause an extremely large t c changes the global behavior 
dramatically. 

The TC model is defined for arbitrary dimension so 
that an extension to three dimensional systems is straight- 
forwardly performed |3ql , and also gravity or moving walls 
can be implemented ]14| , p7[ without loosing its general- 
ity. A future aim is to include the TC model into kinetic 
theories in the style of Haff (|] where it will affect only 
the energy dissipation rate via a correction factor 1 — c e . 
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However, besides excluded volume and dissipation, one 
important property of granular materials, namely friction, 
is missing in its present form. A proper definition of fric- 
tion in the framework of the TC model is in progress and 
has to be tested with systems like particles on an inclined 
plane or sandpiles which are kept at rest by friction. 
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